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1 .   INTRODUCTION 

This  paper  is  concerned  with  the  stability  over  the  finite 
interval  ts[0,T]  of  multistep  formulas  subject  to  varying  step  sizes. 
Although  practical  algorithms  vary  the  step  size  in  order  to  use  as 
large  a  step  as  possible  consistent  with  a  local  control  of  the  estimated 
error,  the  existing  literature — with  the  exceptions  noted  below — deals 
only  with  the  stability  of  fixed  step  formulas. 

The  standard  derivations  of  multistep  formula,  e.g.  Henrici  [3], 
are  based  on  equal  intervals.   There  are  two  common  ways  of  handling 
variable  steps  for  multistep  formula.   One  is  to  interpolate  through  the 
computed  points  to  get  values  at  equally  spaced  points  for  use  in  a 
standard  equal  interval  formula.   The  other  is  to  derive  multistep  formula 
directly  based  on  unequal  intervals .   These  techniques  do  not  have  the  same 
stability  characteristics.   In  1962  Nordsieck  [k]   noted  an  apparent 
instability  in  his  interpolation  version  of  Adams  formula  if  the  step  size 
was  varied  too  frequently.   In  1969  Piowtrowski  [5]  proved  that  a  variable 
step  form  of  the  Adams-Moulton  formula  was  stable  and  convergent.   Brayton, 
et.al.  [l]  and  Gear  [2]  (pp.  1U5-1U6)  point  out  that  the  two  techniques 
for  changing  step  size  are  not  equivalent,  and  Brayton,  et.al.  give  test 
examples  that  show  that  both  can  cause  instability,  but  that  the  technique 
of  step  variation  used  by  Piotrowski  might  be  more  stable. 

The  stability  of  a  method  is  affected  by  three  factors:   the 
underlying  multistep  formula,   the  technique   used  to  handle  variable 
steps,  and  the  scheme   used  to  select  step  sizes — the  formula,    technique, 
and  scheme   for  short.  Method   will  be  used  to  refer  to  the  combination 
of  a  formula  and  a  technique.   This  paper  introduces  the  idea  of  a  step 
selection  scheme,  and  defines  the  stability  and  convergence  of  a  method 


with  respect  to  a  scheme.  A  stability  condition  is  derived.   It  is 
necessary  for  stability,  and,  together  with  consistency,  implies 
convergence  and  stability. 

It  is  easy  to  generate  consistent  formulas.   The  more  difficult 
task  is  to  determine  whether  a  combination  of  formula,  technique,  and 
scheme  is  stable.   The  two  common  techniques  for  step  changing  are 
examined,  and  it  is  proved  that  the  variable  step  technique  is  superior 
to  the  interpolation  technique  for  Adams  formula.   Specifically,  it  is 
proved  that  the  variable  step  technique  used  with  Adams  formula  is 
stable  and  convergent  with  respect  to  any  step  selection  scheme,  while 
the  interpolation  method  can  be  unstable  if  the  steps  are  allowed  to 
change  too  much.   However,  it  is  proved  that  the  interpolation  technique 
with  the  k-step  Adams  formula  is  stable  if  the  step  selection  scheme 
guarantees  that  the  step  is  constant  for  at  least  k-steps  after  each  change 

The  investigation  is  extended  to  all  multistep  formulas,  whether 
explicit,  implicit,  or  predictor-corrector.   It  is  proved  that  if  the 
underlying  formula  is  stable,  then  both  step  changing  techniques  lead  to 
stable  methods  if  the  rate  of  change  of  the  step  size  is  small.   It  is 
also  shown  that  common  error  control  algorithms  lead  to  small  rates  of 
change  in  the  step  size. 


2 .   BACKGROUND 

The  two  step  changing  techniques  will  be  called  interpolation 
and  variable  step,   respectively.   Intuitively,  we  define  them  as  follows. 
Suppose  we  are  using  a  k-step  formula  for  integration  and  suppose  that 
an  approximation  y   .  to  the  solution  is  known  at  k  values  of  the 
independent  variable  t,  say  t   .,  k  >  i  _>  0.   The  interpolation 
technique  for  changing  step  size  to  h  consists  of  the  following  two  steps: 

1.  Interpolate  through  the  known  points  to  get  the  values 

of  the  solution  at  the  points  t   .=t   -ih,k>i>0. 

n-i    n  — 

2.  Use  the  fixed  step  integration  formula  on  these  new 

points  t   .  to  find  the  value  at  t  -,  =  t  +  h. 
n-i  n+1    n 

In  the  variable  step  technique  we  start  with  k  unequally  spaced 
points  t   .  ,  k  >  i  >_  0  and  compute  the  coefficients  of  the  multistep  formula 

Vl  =  "l,n  ^n  +-  •  -+  ak,n  yn-k+l  +  B0,n  \  Cn  (l) 

+  Bl,nh„-lyn+-'-+6k,nhn-^n-k+l 

so  that  it  is  exact  to  the  required  order,  where  h  =  t  _  -  t  .   Since 

n    n+1    n 

the  order  cannot  exceed  2[k/2]+2  if  the  formula  is  to  be  stable  for  fixed 
step  sizes,  additional  restrictions  must  be  imposed.   These  consist  of 
specifying  the  values  of  some  of  the  coefficients  a  and  g.   The  fixed  step 
formula  underlying  the  variable  step  technique  is  said  to  be  the  formula 
that  is  obtained  from  (l)  when  equally  spaced  points  t   .  are  used.   For 
example,  if  we  required  that  a.   =  0,  j  >_  2,  and  that  (l)  have  order  k+1, 
then  Adams-Moulton  is  the  underlying  fixed  step  formula.   If  we  also  ask 
that  g    =0  and  lower  the  order  to  k,  we  get  a  variable  step  method  based 
on  the  Adams -Bashforth  formula. 


The  stability  of  an  integration  formula  subject  to  a  step 
changing  technique  may  depend  on  the  way  in  which  the  steps  are  varied. 
Therefore,  we  must  define  stability  with  respect  to  a  particular  step 
selection  scheme.  First,  therefore,  we  state 

Definition  1       A  step  selection  scheme  is  a  function   9  such  that 

h     =  hQ(t  }h) 
n  n 

where,  for  all  h  >  03   0   <  t   <  T 

1   >  Q(tjh)   1  A  >  0 

NOTE:   As  h  is  reduced  to  zero,  the  maximum  step  size  is  reduced  to  zero. 
The  lower  bound  A  serves  to  prevent  zero  step  sizes  for  any  non- 
zero h;  these  would  not  give  useful  numerical  integration  schemes! 
In  this  paper  we  will  be  concerned  only  with  differential 
equations  y'  =  f(y,t)  for  which  f(y,t)  is  Lipschitz  continuous  in  the 
strip  |y|  <°°,  0  <_  t  <_  T  <  °°.   Small  perturbations  to  the  solutions  of 
such  equations  at  a  given  point  cause  small  perturbations  later  on. 

Stability  means  that  small  numerical  disturbances  at  one  point  t  cause 

n 

small  perturbations   in  the  numerical  solution  at   later  times   for   all 
small  h.      Convergence  means   that   the  numerical  solution  can  be  made 
arbitrarily  accurate  ash+0. 

Definition  2     A  method  is  stable  with  respect  to  a  scheme  6   if  there 

exists  a  constant  M  <  °°   (dependent  on  the  differential  equation  only)  such   that 

\y      -ii        <  M    \y      -  ij 
'^m      vm>  li7n       an] 

for  all  0  <_  t     <   t     <T3   where  y  .  and  y.  are   two  numerical  solutions. 

I  I  III  Is  (s 


Definition  3       A  method  is  convergent  with  respect  to  a  scheme   8  if  the 
computed  solution  y     converges  to 
the  starting  errors   tend  to  zero. 


computed  solution  y     converges  to  y(t  )  for  any  0  <_  t     <_T  as  h   ■*■  0  and 


We  'will  only  discuss  a  single  differential  equation.   The 
extension  to  systems  of  equations  is  straightforward,  although  notationally 
complex.  We  will  also  assume  that  f(y,t)  and  y(t)  are  differentiate 
as.  often  as  necessary  to  save  stating  all  the  conditions  each  time.   The 
results  are  valid  if  f  is  only  Lipschitz  continuous  ,  and  can  he  obtained 
"by  replacing  f  =  3f/9y  by  the  Lipschitz  constant  L  whenever  f  appears. 


3.   NOTATIONAL  DEVELOPMENT 


By  developing  a  uniform  notation  for  all  the  formulas  and  step 

changing  techniques  to  be  discussed,  we  will  be  able  to  give  general 

proofs  in  Section  h   that  can  be  applied  to  all  cases .   Some  of  this 

notation  has  been  developed  in  [2],   It  is  summarized  briefly  here  and 

suitably  extended. 

Consider  the  integration  formula  (l)  with  3„   =0.   It  is  an 

0,n 

explicit  formula  which  can  be  used  directly.   It  can  also  be  used  as  a 
predictor  for  a  corrector  formula.   We  will  write  the  corrector  formula  as 


*M]    =  ^n  ^n  +-"+  «k,n  ^n-k+l  +  gS,n  hn  f(Ci>  W 

+  (B*   h  .  y'  +...+  g*   h  ,  y«  ... 

l,n  n-1  n        k,n  n-k  n-k+1 


(2) 


.(0) 


A  predictor-corrector  formula  is  one  in  which  a  predicted  value  y    is 


computed  by  the  explicit  equation 

.(0) 


y'n=ot,    y  + . . .+  a    y  ,  , 
Jn+1    l,n  °n        k  ,n  ^n-k+l 


(3) 


+  3n   h    y'  +...+  e    h    y'  ^ 
l,n  n-1  n        k,n  n-k  n-k+1 


and  the  equation  (2)  is  used  to  compute  a  fixed  finite  number  of  corrected 

(l)       (M)  (m) 

values  y  ,.,...,  y   '.   The  approximation  y  ,.  is  defined  as  y  , '  and  y',. 
n+1       n+1  n+1  n+1      n+1 

is  set  to  either  f (y    )  or  f(y     )  depending  on  whether  or  not  a  final 
function  evaluation  is  needed. 

An  implicit  formula  is  one  in  which  equation  (2)  is  solved  for 
y      =  y    .   One  way  of  doing  this  is  to  iterate  a  predictor-corrector 
formula  until  successive  iterates  are  approximately  equal.   This  converges 
for  h  sufficiently  small  but  may  diverge  for  large  h.   Other  methods,  such 
as  Newton-Raphson,  can  also  be  used  to  solve  (2).   The  stability  of  a 
method  based  on  an  implicit  formula  is  independent  of  the  process  used  to 
solve  equation  (2),  so  we  do  not  need  to  specify  that  process. 


During  the  computation,  a  number  of  previous  values  of  the 


function  y   .  and  the  scaled  derivatives  h   .  _  y1  .  must  be  saved, 
'n-i  n-i-1  n-i 

Let  us  collect  these  together  in  a  column  vector  y  defined  "by 


^n  =  [V  yn-l: 


»  K  ,  y 1 ,  K  o  y '  -. » •  •  • »  *   t.  y '  ^  ] 


T 


n-k+1'  n-1  Jn*  n-2  Jn-1 


n-k  ' n-k+1 - 


for  a  k-step  formula,  where  T  is  the  transpose  operator.   (if  some  of 
these  values  are  not  needed,  we  can  drop  them  from  the  vector.   For 


example,  in  Adams  formula,  y  would  be  [y  ,  h  _  y ',...,  h  .  y'  .  ,. 

*-n  n   n-1  n       n-k  n-k+1 

A  numerical  integration  method  is  a  method  for  computing 


]T.) 


y  from  y    .      We  will  first   describe   it   for  the  variable  step  technique, 

Let  y\,i  be   defined  as 


r    (m) 


y:i?}>  K  n  y:.-.-» 


n+1'   *&""*  °n-k+2'      n  Jn+1    '      n-1  Jn 

where  h     y'*?     is  h     f(y(™,       ,  t    ..  )   for  m  >  1, 

n     n+1  n      wn+l  n+1  — 


y'  1 

n-k+1  <yn-k+2J 


hn  yn+l  Yl,n  yn  +  Y2,n  yn-l  +-'-+  \>n  yn-k+l 


+   6  h  v'    +   5  h  y'         +...+  6  h  y' 

l,n     n-1  ^n  2,n     n-2  °n-l  k,n     n-k  ^ n-k+1 

when  m  =  0,  and  the  coefficients  y   and  6  are  given  by 

Y-    =  (a.    -  a*   )/6* 
i  ,n     i  ,n    i  ,n   0,n 


and 


6.   =  (6.    -  (3*   )/6* 

i,n     l  ,n    i,n   0  ,n 


From  these  definitions  we  can  write 


til  -  *n  *, 


(h) 


where  A  = 
n 


— 

— 

l,n 

2,n 

•      •      • 

k-l,n 

k,n   j    l,n 

B2,n 

•    •    • 

ek-l,n 

6k,n 

1 

0 

0 

0 

0 

1 

0 

0 

0 

• 

•• 

• 

0 

0 

1 

0 

l,n 

Y2,n 

9    ■     • 

Yk-l,n 

Y            !  5 

k,n   j    l,n 

62,n 

6k-lsn 

k,n 

1 

0 

0 

0 

0 

0 

1 

B     •      • 

0 

0 

• 

•    • 

•    •    • 

• 

. 

i     o 

0 

1 

0 

J 

and 


(m+l)    (m)     r       f(y(m),  t   )  -  h  y'U)] 
^n+1     ^n+1   V  n  K7n+V     n+r    n  yn+l  J 


(5) 


where 


*   =  t3*   ,  0,...,  0,  1,  0,...,  0 


(6) 


Equation  {h)   represents  the  predictor — the  evaluation  of 

y  ,,  and  h  y'  Ln  as  linear  combinations  of  the  known  information  plus 
"n+1      n  J    n+1  ^ 

the  shifting  over  of  the  saved  information  and  the  discarding  of  y 

and  y'    .   Equation  (5)  represents  the  corrector — the  evaluation  of 

s-t    (m)  a.        \        a  j.-,  n   ^  •  j>     (m.)    ^  •,    ,(m)   _ 

£\Y   xD  "t  ,  -,  )  and  the  updating  of  y   '  and  h  y ',.  .  For  convenience, 
n+1   n+1  n+1      n  n+1 

we  define 


/  (m)  x   .   _/  (m)   ,    >.   .    ,  (m) 


V^n+1 


n+1'  n+1' 


n  "n+1 


(T) 


Then  equation  (5)  can  "be  written 

(m+l)    (a)+l   F  (  (m)}  (8) 

^+1     *n+l   -n  nv,Ln+l;  v  ; 

If  a  predictor  only  formula  is  used,  we  still  formally  must  apply  one 

step  of  (8)  with  H      =  £_  =  [0,. .  .,  Q,  1,  0,,..,0]  .   (in  this  case  the 

value  of  a*  and  3*  are  unimportant.)   If  a  final  function  evaluation  is 

used,  we  apply  (8)  M  times  using  the  &   given  in  (6),  and  then  apply 

p 
one  further  step  of  (8)  using  &   =  l_   . 

Now  let  us  consider  a  fixed  step  method  in  preparation  for 

formalizing  the  interpolation  technique.   In  this  case,  the  a,  3,  y,  and 

6  in  A  and  I      are  independent  of  n.   Suppose  we  have  just  completed  the 

step  from  t   .  to  t  using  step  h   . .   The  values  saved  in  y  are  the 
n-1     n  n-1  •4i 

computed  approximations  at  the  points  t  -ih    ,  0  <_  i  <  k.   Before 

stepping  to  t  , ,  =  t  +  h  we  must  interpolate  to  the  points 
**  n+1    n    n  r 

t   .  =  t  -  ih  .   This  is  a  linear  process  and  can  he  represented  by 
n-i    n     n 

first  premultiplying  y  by  an  interpolating  matrix  C  .   Consequently, 
we  write  the  predictor  step  for  the  interpolation  process  as 

C\  -  A   Cn  *n  (9) 

The  corrector  step  (8)  is  unchanged  except  that  I      =  l_  is 
independent  of  n. 

Both  for  theoretical  development  and  computational  ease,  it  is 
convenient  to  perform  a  transformation  based  on  Nordsieck's  [k]   form  of 
Adams  formula  for  the  interpolation  technique.   This  consists  of  computing 
and  saving  Q  y  rather  than  y  where  Q  is  a  nonsingular  matrix.   The 
particular  form  of  Q  used  is  obtained  as  follows:   The  set  of  values  of 
y  _.  and  y'_.  0  <_  i   <   k  determine  a  unique  (2k-l)-th  degree  polynomial. 


10 


Represent  that  polynomial  instead  by  its  value  at  t  and  the  values  of 
its  first  2k-l  scaled  derivatives  hJ   y   /j  !  ,  1  <_  j  <  2k.   (This 
leads  to  what  is  called  a  2k -value  method  in  Gear  [2].)   If  y  is  the 


set  of  saved  values  of  y  and  hy'  "based  on  step  h 


n-11 


and  a  is  the  set 
— n 


(10) 
(11) 


of  values  of  [y  ,  h  .  y',  h2  _  y"/2,...,  h2*"1  y(2k-l)/(2k-l) ! ]T 
n   n-1  n   n-1  n        n-1   n 

representing  the  same  polynomial  of  degree  2k-l,  then  Q  is  defined 
"by  a  =  Q  y  •   Q  is  independent  of  y  and  h.    .   If  we  apply  this 
transformation  to  (9)  and  (8)  we  get 

Q  4°|  =  (QAQ"1)(QCnQ~1)Q  y^ 

To  simplify  the  presentation  we  will  write  y  rather  than  a  ,  A  rather 

than  GAQ   ,  C  rather  than  QC  Q   ,  I      rather  than  Q£  ,  and  F(c)  rather 
n  n     — n  -n        — 

than  F(Q  c_)  when  we  are  talking  about  the  transformed  scheme. 
Consequently,  equations  (10)  and  (ll)  are  identical  to  (9)  and  (8). 
The  advantage  of  making  the  transformation  is  that 


pv_1 
Cn  =  diag[l,  hn/hn_l5...,  (\/\_±)  1  (12) 

so  that  the  interpolation  process  is  simplified.   This  was  Nordsieck's 

motivation.   Finally,  we  note  that  (9)  is  identical  to  (k)    if  A  is 

n 

interpreted  as  A  C   for  the  interpolation  technique.   Consequently, 
equations  (k)   and  (8)  represent  the  predictor  only  formula  or  the 
prediction-correction  formula  in  either  step  changing  technique. 

The  corrector  only  formula  can  sometimes  be  implemented  by 
iterating  the  corrector  to  convergence  and  treated  theoretically  by 
setting  y    =  lim  yV  ,  •  However,  this  process  may  not  converge  (as 
in  stiff  equations),  in  which  case  the  corrector  must  be  solved  by  other 


11 


means.     We  note  that  the  corrector  equation   (8)   always  adds   on  a 
multiple  of    &     to  the  current  value.      Consequently,   the  final  value 
has   the  form 

Vl  ■  J&i  +  "4  =  An  *n  *  «*»  (13) 

where    W  is  a  scalar.      The  value  of    go  is  such  that 

nzlH  *  co^)  =  o  (no 

(Although  we  have  used  predictor  coefficients  in  A  ,  some  straightforward 
arithmetic  will  show  that  the  final  value  of  y  is  independent  of  these 
coefficients  when  the  corrector  equation  is  solved  by  (13)  and  (lU).) 

We  must  now  develop  the  error  propagation  formula  for  equations 
(k)    and  either  (8)  or  (13)  and  (lU).   The  global  error  n  is  defined  by 

*n  "  ^  -  £(ta)  (15) 

where  y_(t  )  is  the  value  of  the  components  of  the  vector  y  as 
evaluated  along  the  desired  solution  of  the  differential  equation. 

The  truncation  error  d  is  the  error  introduced  by  one  step  of  the 

-n 

method,  starting  from  the  solution  of  the  differential  equation.   Thus, 
d  is  defined  by  the  equations 


either 


&?  -  tit  ♦  4  FnO  - 


^n+1       ^ti+I 


if  M  corrector   iterations   are  used,   or 


2^4.t    =  liVi    +  wi     =  A     y_(t    )   +  Sd      ,  (16) 

*-n+l       *-n+l  — n  n  *-    n  —a 


12 


where  F(jf  _  )  =  0  if  the  corrector  is  solved  exactly,  and 

4  =  4+1  -  2L(tn+1) 

From  the  definition  of  e  , ,  and  d 

-n+1  -11 

^n+1  =  4+1  "  1.(^+1 ) 

"  4+1   "  4+1  +  4+1  "  ^    n+r 

=  2,,+!  "  4+1  +  4  (1T) 

(Note  that  if  round-off  errors  are  to  "be  considered,  they  can  be  included 
in  d  at  this  point,  so  that  the  computation  of  y  _  can  be  considered 
to  be  exact.)   If  M  corrector  iterations  are  used,  (17)  can  be  written  as 

-n+1   ^+1   *ia+l   -n 

=ei)+^vei))-ei)-4vaii,>^ 

3F(M) 
by  the  mean  value  theorem.   - —    is  the  row  vector  of  partial  derivatives 

of  F  with  respect  to  the  components  of  y_  evaluated  somewhere  in  the 

interval  (y     ,  f  ).   From  equation  (7)  we  see  that 

|^M  -  [h  f<M),  o,...,  0,  -1,  0,...,  0] 
dy_       n  y 

where  rM'  is  9f/3y  evaluated  in  the  interval  (y^~   ,  y  +~   ) .   The 

process  leading  to  (l8)  can  be  repeated  to  get 

,(M)  _(D 


^n+1 


-n  3y_          -n  3y_     "-n+l  *-n+l   -n 

aF(M)  .  (1) 

-  [I  +£   I1   ]...[!  +£   I1   1  A  e  +  d               (19) 

-n  3y_           — n  3y_      n  -n  n 


We  write  this  in  the  form 
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e  Ln  =  S   e  +  d  (20) 

-n+1    n  -n   -n 


T 
Let  e.  be  the  unit  row  vector  consisting  of  a  one  in  the  i-th 
—i 

T 
position  and  zeros  elsewhere.   Thus,  e_  is  [l,  0,...,  0].   Let 

subscript  d  represent  the  component  corresponding  to  the  h    jl 

entry  in  y  .   (d  is  k+1  for  the  untransformed  method,  2  for  the 

transformed  method.)   Then 

fl(m)  .  h  f<»>  J  .  £  '    (21) 

9v_       n  y   — 1   -d 

Therefore,    from   (19),    (20),    and   (21) 

's     =      n      [(I   -I      e?)   +  h     f(m)   I      e?]A  (22) 

n  .  -n-d  ny-n— In 

m=l 

T  M 
=  [(I  -  I      e, )   +  I      {Products  of  M  terms  of  form 
—n-d 

h   f(m)  I      e%r  (I  -  I      eT) 
n  y   -n  — 1        -n  — d 

with   at  least  one  of  the   former} ]A 

n 


We  define 


Then  we  can  write  (22)  as 


/        T.M  . 

S  =  (I  -  I      e,uA  23 

n       -n  — d   n 


S  =  S  +  h   S  (21+) 

n    n    n  n 


where  S  is  a  matrix  whose  elements  are  polynomials  in  h  ,  f    and 
n  n   y 

entries  from  A  and  %    . 
n     — n 


Ik 


Lemma  1   If  a,  g,  a*  and  3*  are  uniformly  "bounded  or  if  A  and  I 
n     -n 

are  uniformly  "bounded  so  is  j  |S  |  j  for  all  |h|  <_  u. 


Proof     The  elements  of  S  are  polynomials  m  h  whose  coefficients 
_ n  n 

are  polynomials  in  f    and  the  entries  of  A  and  I     .   If    I  is  bounded 

y  n     -n     y   ' 

by  the  Lipschitz  constant,  while  entries  in  A  and  I     ,  except  possibly 

n     — n 

for  the  d-th  row  of  A  ,  are  bounded  if  a,  3,  a*  and  3*  are  bounded.   The 

n 

d-th  row  of  A  contains  elements  of  the  form  (a-a*)/3*   ,  and  could  be 
n  0  ,n 

unbounded  if  3*   can  approach  0.   However,  from  (22)  and  (2U),  the  only 

way  the  d-th  row  of  A  enters  S  is  when  A  is  premultiplied  by 
*  n         n         n 

T 
(I  -  i      e ., )  which  zeros  the  d-th  row  and  changes  the  first  row  from 
— n  -d 

a  to  a*.   Consequently,  the  coefficients  of  h   in  S   are  bounded,  so  for 

*  n     n 

any  u  >  0,   Is  ;   is  bounded  for  all  0<h<u.   If  A  and  I      are 
j  •>    i  i  ni  i  _   _        n     -n 

uniformly  bounded,  the  result  follows  similarly. 
Q.E.D. 


The  behavior  of  the  error  e   is  determined  by  equation  (20). 

-n 

When  I  Is    is  bounded,  we  can  use  equation  (2U)  to  relate  the  stability 
n 

of  the  method  to  the  properties  of  S  .   The  following  condition  will  be 


shown  to  be  necessary  for  a  method  to  the  stable: 
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Definition  4      A  method  satisfies   the  stability  condition  with  respect 
to  a  step  selection  scheme  if  there  exists  a  C  <  <*>3   such  that 

n         m-1     m-2  n 

satisfies 

|  1^1  |  <_C  (26) 


n' 


for  all  0  <   t     <  t     <  T. 
J  —    n        m  — 


If  an  implicit  method  is  used,  the  corrector  is  solved 

exactly,  and  we  must  derive  equations  (20)  and  {2k)   by  a  different 

technique.   Substitute  (13)  and  (l6)  in  (17)  to  get 

e   n  =  A   e  +  (w-fiS)  i    +  d  (27 

-n+1    n  — n         —   -n 

(The   subscript  n  has  been  omitted  from  l_  for   clarity  later.)     We  know 

0  =  F(vii+1)    -  F(2n+1) 

9v_  ^+1       *n+l' 
3F 


3v_  <An  ^  +  («-a)  D 


Hence , 


(u-fi)  =   -C|Jl>        #A     e    } 

3v_  —  3v_    n  -n 

Note   that  both  quantities   enclosed   in  braces   are   scalars.      Substituting 
in   (27)  we   get 

S„+i    =Ae      "   ^^>       l    <■¥■  A     £    >   +  d 
-n+1  nn  9y_—         —     3v_    n  -n  -n 

°%_  ~         - "  ov_       n  -n       -n 

=  S     e      +   d 
n  — n       — n 
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where 


n         9y_  -    -  3y_   n 

i-  [I  -  "'{L  h  f  -  a}"1  £(h  f  e^  -  e^)]  A 
1  n  y    d    —  n  y  — 1   -d/Jn 


Here,  I  .    is  the  i-th  component  of  l_. 

If  we  choose  u  such  that  for  0<h  <  u,  h  JL   f  -  £   is  hounded  away 

n  —  n     1     y  d  J 

from  zero    {l      is   one   in  the  methods   discussed),   then  we  can  write 
|hn'l   fy  "V"laS   ^f   {1  +llhn  fyAd  +°(hn>>-      He^, 

sn=  [I-^Ajg]  An  +  hnsn  '        (28) 

=  S     +  h     S 
n  n     n 

The  form  of  S  is  the  same  as  the  previous  form  with  M  =  1  and  I    scaled 
n  — 

so  that  Jl   =  1,  as  it  always  is  in  methods  discussed  in  this  paper. 
We  also  need  to  extend  Lemma  1  to  this  case. 


Lemma  la  ||  S  ||  defined  in  (28)  is  hounded  if  a,  6,  a*  and  B* 

or  A  and  I      are  uniformly  bounded, 
n     -n 


Proof     From  (28)  and  the  preceding  equation 

S     =   -Unh   f  -I  A'1  f     I    e.    A 
n  lnyd  y In 

+    [U.h   f  -4  J"1  +   C"1]^    eT  A  /h 
lnyd  M dnn 

=   -Unh   f  -2,  J-1   f    I  e?  A 
lnyd  y In 

+  nT^H-h   f  4)"1   f     a   s,      e^  A 
d        lnyd  y—     1-dn 


IT 


All  "but  possibly  the  d-th  row  of  A     are  bounded.      The  first  term  does 

not   involve  the  d-th  row  of  A   ,   and  is,   therefore,  bounded  for 

n 

0  <  h  <  u.   The  second  term  involves  the  d-th  row  of  A  ,  but  scales  it  by 
—   —  n 

£ .,  =  3*   .   Hence,  it  is  also  bounded  for  0  <  h  <  u.   Hence,   Is    is 
1    0,n  —   —  '  '  n1  ' 

bounded.   If  A  and  I      are  bounded,  the  result  follows  directly, 
n     — n 


Q.E.D. 
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h.      STABILITY  AND  CONVERGENCE 

Theorem  1  If  a  method  is  stable  with  respect  to  a  step  selection 

scheme   6j  it  satisfies  the  stability  condition   (26)   with  respect  to 
that  scheme. 

Proof     Consider  the  differential  equation  y'  =  0  for  which  S  , 
n 

defined  "by  (2k),    is  zero.   Evidently,  the  difference  between  two 
numerical  solutions  satisfies 

^n+1  ~  -n+1  ~  Sn^n  "  Za) 

i  i  Am  i  i 
If  the   S  II  are  not  uniformly  bounded,  there  exist  t  and  t   and 

1  '  n1  '  m      n 

bounded  y  -  y  such  that 
— n   — n 

ym   ~  Jm  =   Sm(y  -  y  ) 

is  arbitrarily  large,  implying  that  the  method  is  not  stable  with 
respect  to  0.   Hence  stability  implies  the  stability  condition. 
Q.E.D. 

The  convergence  theorem  depends  on  the  following  lemma. 

Lemma  2   If,  for  the  difference  equation, 

x^=Sx+hSx+hA  (29) 

—n+1    n  -n    n  n  -n    n  — n 

there  exist  constants  k  ,  k  ,  and  k  such  that 
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Is  |  |  <  k.      ¥  n 
n   —  1 


A    <  v       ¥  n, 
—n   —  2 


and  |  |sm|  |  <  k.      ¥  m  >  n 

n   —  0         — 


where 


sn  =  I 

n 


and  S   is  defined  by  equation  (25) a 
n 

then,  if  k   >  0 


or  if  k  =0 


I  II   rv  II  II   k2i  Vl(W   k2         ,,M 

1^11  1^  M^ll  +-]  e  --  (30) 


1^1 1  i^o  H*qII  +v2(Vto)  (31) 


Proof     From  (29)  we  have 


N-l 

x=  s'  x^  +  E   s  4.1  h  S  x  +  A 

— N    0  -H)     „   n+1  n  n  — n   -n 
n=0 


Therefore , 

112,11   ik0    I  I*,  I  I    ♦     I     Vn(kl   1 1^1 1   +  k2>  (32 

n=0 

Consider  the  case  k  >  0  first.   We  show  (30)  by  induction.   It  is 
evidently  true  for  N  =  0  since  k  >_  |s  |   =  1.   Assuming  its  validity 


for  n  <  N,  we  substitute  in  (32)  to  get 
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N-l  k„       knkn(t  -tA) 

I  I  xl  I    <  k      11x11+     E      kkh    fk      llxll+— I*-01     n      ° 
I  l%ll    ±kQ    |  IXqI  I    +     L     KQK  nnLKQ    HXqI  I    +  k   J   e 

n=0  1 

k0  N-l  k  k    (t   -tj  k„ 

n       II       II  2irn    j.     v      i    i    v.  0   1     n     0    t  2 

=    [kQ    MxqII    +— ][l  +     E     Wn  e  ]   -  — 

1  n=0  1 

k0  N-l       knk,h  k_kn(t  -tj  k0 

^    n       II       ii     .      2  T  r  _,  _,/01n_x01vn0/1  2 

1  n=0  1 

This  completes  the  inductive  proof  of  (30).   If  k  =  0,  equation  (31) 

follows  from  (32)  directly. 

Q.E.D. 


The  convergence  theorem  is  now  straightforward. 


Theorem  2  If  a  method  satisfies   the  stability  condition  with  respect 

to  a  step  selection  scheme   0,  if   a,  a*.  3  and   6*  or  A     and  I     are  uniformly 

bounded,   and  if  the  truncation  error  d    satisfies 

3  J  — n 

I  Id   I  I    <  h     e(h)        V  n 
i \_n\  i    _    n 

where  e(h)  +  0  as  h  -*■  03    then  the  method  is  stable  and  converges  with 
respect  to  <d3   and  there  exist  constants  k~  and  k.  such  that  the  global 
error  e_^  satisfies 

llf^ll   £k3    ll^ll    +  k4  e(h)  for  0  <_tN<_T 


Notes:  1.      The   condition  on  d     is   a  consistency  condition.      If  the 
-n 

method  has    order  r,   then  e(h)    =  0(h    ).      However,    a  more  useful  way   of 

looking  at  this  result    is   that   if  we   control  the  step   size  h     so  that   the 

n 
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truncation  error  is  "bounded  by  eh  ,  the  global  error  is  bounded  by 

k   i  |e  |  j  +  ki  e,  so  that  changes  in  the  truncation  error  control  are 

proportionately  reflected  in  the  global  error  bound. 

2.   The  condition  that  a,  a*,  3  and  $*  or  A  and  I     be  uniformly 

n     -n 

bounded  is  also  a  form  of  consistency.   It  rules  out  the  following  type  of  case 
The  explicit  formula 


1 


n+1   h  „-h  . 
n-2  n-1 


(h  +h   J2-h2  .  (h  +h   J2-h2  . 

n  n-1    n-2  n  n-1    n-1 

h   +h  _   ~  yn  h  0+h  .   "  Jn-2 

n-2  n-1  n-2  n-1 


h  (h  +h   +h   _) 
n  n  n-1  n-2 


hn-2       ****■ 

is  second  order  provided  that  h   _  ^  h   .  .   However,  as  h   _  ■*■   h   .  ,  the 

n-2    n-1  n-2    n-1 

coefficients  and  the  error  term  blow  up.   If  we  used  this  formula  with 

a  step  selection  scheme  which  kept  h   _/h   .  bounded  away  from  one,  we 

n-2     n-1 

would  get   convergence   if  the  method  was    stable,   but   if  the  step  selection 

scheme  let  h     . /h     „  •*  1  as  h  -*-  0 ,  we  may  not   get   convergence. 
n-1     n-2 


Proof  From   (20)    and   (2U)    the   error   equation   for   a  method  is   given  by 

e    ...    =  S     g     +  d  (33) 

-n+1         n  -n       — n 

■  S     e      +h     S      e     +d 
n  — n  n     n  — n       — n 

Consistency  and  Lemma  1  imply  that   |s  |  <^  k.  for  some  k  <  °°  provided 

that  h   <_  u.   Lemma  2  can  be  applied  to  (33)  with  k  =  e(h)  since  the 

stability  hypothesis  guarantees  the  existence  of  a  constant 

k^  such  that   I S  II  <  k^ . 
0  '  '  n '  '  —  0 
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Therefore, 

1  k3  I  1%  I  I  +  \   e(h)  v  °  1  tN  1  T 
Hence,    e^ J  |  -*■   0  if   |  e^  |   and  e(h)  -*■   0  and  we  have  convergence. 

— Vi  — 0 

The  difference  between  two  numerical  solutions  satisfies 

^n+1  "  4+1  =  Sn(yn  "  V 
The  same  analysis  (with  e(h)  =  0)  yields 

Therefore,  the  method  is  stable. 
Q.E.D. 

The  difficult  problem  left  is  to  verify  the  stability  condition 
for  various  methods  and  step  selection  schemes.   First,  we  will  deal  with 
the  Adams  methods  in  Theorems  3,  h,    5,  and  6.   Theorem  5  is  an  extension 
of  the  result  of  Piotrowski  [5]. 


Theorem  3  A     and  i      are  bounded  for  the  interpolation  technique 

using  a  step  selection  scheme   0. 


Proof     For  the  interpolation  technique  I      =  i_   is  independent  of  n,  and 

is,  therefore,  bounded.   A  =  A  C  where  C  ,  given  by  equation  (12),  is 

n      n        n 

— 2k+l 
bounded  by  A     "  from  definition  1.   Therefore,  A  is  bounded. 

n 

Q.E.D. 
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Theorem  4  a,   a*,   3  and  3*  are  bounded  for  the  variable  step  Adams 

method  using  a  step  selection  scheme  9. 

Proof  The  variable   step  Adams  method  has  the   form 

y  =  v     +  R*       h     v'        +...+  3*       h         v1  (35) 

Jn+1       *n        P0,n     n  ^n+l  Pk,n     n-k  ^n-k+l  u?i 

where  the  3*   are  determined  from  the  requirement  that  the  method 
i,n 

have  order  k+1.   (For  an  explicit  method,  3*   =0  and  the  order  is  k.) 

0  ,n 

A  necessary  and  sufficient  condition  for  the  order  to  he  k+1  is  that 

(35)  he  exact  for  the  functions  y(t)  =  (t  -  t    )  ,  0  <_  r  <_  k+1.   This 

leads  to  a  system  of  k+1  linear  equations  for  the  3*   ,  namely 

i  ,n 

k 

Z  3*   h  ,  r(t  ..  ,  -  t  _,,)  (36) 

.  _   i,n  n-i    n+l-i    n+1 
i=0 

=  -(-hn)r,  r  =  1,  2,...,  k+1 
Define 

Vl-i  "  tn+l 
co .  = 


i         h 


h  +h  .+...+  h  ..  . 
n    n-1        n+l-i 


From  definition  1 


and 


A<_|oo.|<i,    i  >_  1  (37) 

0  =  Uq 


oj.-u.I  >  A      if  i  4   J  (38) 


Divide  (36)  by  rh  to  get 


,  r 

Z  3*.  (w.-ok..  )uk   =  -  —                (39) 

.  _   i,n  1  l+l  1       r 
i=0 


2k 


The  right  hand  side  of  (39)  and  the  coefficients  of  s*   on  the  left  hand 

i,n 

side  are  bounded  uniformly  in  n  by  (37)  •   Ike  determinant  D  of  the 
coefficients  is  a  simple  multiple  of  the  Vandermonde  determinant,  namely 


D  =  [   n    (w.-oj.)][  n   (w.-o).  ._)] 
k>i>j>p   X   J    i=0   ^     1+1 

By  (38)  this  is  hounded  away  from  zero,  hence  the  solution  of  (39)  is 

bounded  uniformly  in  n . 

Q.E.D. 

Theorem  5  The  variable  step  Adams -Bash  forth    (AB)3  Adams -Bash forth- 

Moulton  predictor-corrector   (ABM)  3   and  Adams -Moulton   (AM)   methods  are 
stable  and  convergent  for  any  step  selection  scheme. 


{ho) 


1 

Proof  The  vector  of  saved  values   y     is    [y    ,  hy',...,   hy '       , . ]      so 
*-n  n  n  n-k+1 


the  matrix  A     is    given  by 
n 


l,n       2,n  k,n 


0     6n  6_  ...      6, 

l,n        2,n  k,n 


while  2,      is    [0,    1,   0,.  ..,   0]      for  AB  and    [3*      ,    1,   0,. 
— ti  u  ,n 

and  AM. 


,   0]      for  ABM 


Hence,  from  equation  (23), 


S   =  (I  -&      eT)M  A 
n        -n  -2    n 
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1 

nl,n 

n2,n 

•  •  • 

\,n 

0 

0 

0 

0 

0 

1 

0 

0 

Ui) 


where 


i,n 


3.    for  AB 
i,n 


i,n 


for  ABM  and  AM 


From  Theorem  h,   n    are  uniformly  bounded,  therefore,  S   is 
i,n  n 


uniformly  bounded.   Hence  the   |S 


n+mi 


are  uniformly  bounded  for 


0  <  m  <  k.   From  the  structure  of  S  ,  we  see  that  for  m  >  k 
—   —  n 


*n+m  __  ~n+k 
n      n 


1  x. 


\ 


(U2) 


Hence  the  I  Is      are  uniformly  bounded  for  all  m  >  n  and  the  method 

satisfies  the  stability  condition.   From  Theorem  k,  A  and  I      are 

'    n  -n 

uniformly  bounded.      Finally,  we   can  expand  y(t      .  )   and  hy'    .    n(t      .) 

n-i  n-i-1     n-i 

by   Taylor's   series   about   t  with   a  remainder   involving  h  y        (?)• 

If  the  order   is   greater  than  r-2 ,   then  all  but  the  remainder  terms 

r— 1  ( r ) i 

will  vanish,   and  we   can  bound       d         by  h    (h  C     Max  y  ). 

1  '-n1  '     n       r     '"    ' 

Thus,  the  three  conditions  of  Theorem  2  are  satisfied. 
Q.E.D. 
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Theorem  6  If  the  step  selection  scheme  is  such  that  at  least  k 

steps  of  constant  size  are  taken  between  step  changes,    then  the  k-step 
interpolation  AB3  ABM3   and  AM  methods  are  stable  and  convergent. 

Proof  Once  again,  we   show  that  the  hypothesis  of  Theorem  2  are 

satisfied.      We  will  use  the  Nordsieck  vector  to  express  the  method  as 

a  k+1  value  method.      Then  the  matrix  S      is   given  by 

n 

§  -  [I  -  A  elf   A  C  (1*3) 

n        z.  n 

=  R  C 
n 

The  matrix  R  has  the  important  properties  that  the  simple  eigenvalue  one 

corresponds  to  the  right  eigenvector  e_  ,  and  that  all  other  eigenvalues 

are  zero  (Gear  [2],  pp.  138-lUo).   Consequently,  R  x  =  y  e_  where 

|y  |/||x||  is  uniformly  bounded  for  all  x  (Tu  [6],  p.  20).   The 

eigenvalue  one  of  C   corresponds  to  the  eigenvector  e_  .   When 

h  =  h  n  ,  C  =1.   From  the  hypothesis  of  the  theorem,  it  is  not 
n    n-1   n  ' 

possible  for  C.  ^  I  ^  C.  if  |i-j|  <  k.   Consequently,  the  product 

Sm  =  RC   nRC   .  ...EC 
n      m-1    m-2        n 

either  contains  less  than  2k+2  matrices,  or  can  be  reduced  to 

Sm  =  R  C   ,  .  .  .  R  C  Rk  C  ,  R^-k_n 
n      m-1        q     q-k 

where 

0  <_  q-k-n  <  k 

Since  the  ratio  of  adjacent  steps   is  bounded  by  A      ,      |C 

ii  mi  i 

is  uniformly  bounded.   Hence,  ||C    R        <_  k  uniformly.   Therefore, 

if  S   contains  2k+2  or  more  matrices 
n 
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Smz=RC     n  ...RC     RkC     .    R1-k-n 
n  —             m-1  q  q-k  n  ± 

=  R  C     n  ...EC     Rk  x 
m-1  q         — 


where      |  x |  |    <_  k^      |  zj  |  . 
Hence, 


m 

n  —  "m-1    ' '  "    il  "q  '"x  — i 


S„    z  =  R  C_  .,    . .  .   R  C_   y„  e 


=   yx% 


so  that 


Smz||    =    Iv 


n  — ' 


x    '  '— ' 


1kg      |x||  (as   y,y||x||    is  bounded] 

lk5  k6    I  111  I 


Therefore,  lsm||  <  k  k, 


n 


^   ~6 


If  S  contains  less  than  2k+2  matrices,   S    is  bounded  by  the  finite 

n  '    n1  ' 

product  of  the  norms  of  each  matrix,   hence  the  method   satisfies   the 

stability   condition. 

Theorem  3  shows   that  A     and  I      are  bounded,      d     is  bounded 

n  -n  -n 

by  the  same  argument  as  used  in  Theorem  5.   Hence  Theorem  2  proves 

convergence  and  stability. 

Q.E.D. 

Tu  [6]  points  out  that  the  interpolation  form  of  the  Adams 
formula  is  not  necessarily  stable  under  any  step  selection  scheme.  He 
shows  that  for  the  three  step  ABM  formula  with 

h2i  =  Wh2i+1  =  h2i+2  =  ^21+3  =  * '  ' 
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the  matrix  S      has  an  eigenvalue  of  [(co-l)  /Uoo]  .   For  sufficiently 
large  w  this  is  unbounded.   Examples  of  instability  when  oj  =  10 
are  given  in  that  thesis  along  with  examples  showing  that  the  variable 
step  technique  is  stable  for  the  same  step  selection  scheme.   This  case 
points  out  that  the  variable  step  technique  is  more  stable  in  some  cases. 

Another  example  in  Tu's  report  substantiates  Brayton,  et.al.'s 
claim  that  the  variable  step  technique  is  more  stable.   The  problem 
y'  =  -y  with  h   =  .05  and  h     =  .005  was  integrated  using  a  three 
step  backward  differentiation  formula.   It  was  apparently  stable  for  the 
variable  step  technique  and  apparently  unstable  for  the  interpolation 
technique.   It  is  hypothesized  that  the  variable  step  method  can  also 
cause  instabilities,  but  no  examples  have  yet  been  found. 

In  practice,  we  do  not  vary  the  step  in  an  extreme  way.  When 
the  step  changes  are  controlled  to  be  small,  both  techniques  are  stable 
as  shown  below. 

Theorem  7  If  a  method  satisfies  the  stability  condition  for  fixed 

stepSj    it  also  does   so  with  respect  to  a  step  selection  scheme   that 
produces  step  size  changes  small  in  the  sense   that 

hn+l 

=  1  +  0(h)   as  h  +  0  (44) 


h 
n 


Note 


If  6(t,   h)    =   6(t)    and  0    is   a  boundedly  differentiable   function 


of  t,   then  the  hypotheses   of  the  theorem  are   satisfied  since 


h    ._         6(t    .,)  6'(£    ) 

n+1  n+1  _  ,    .   .  n 

"h~ =    e(t  )  ~1  +  \  "em 

n  n  n 


=  1  +  h0'(C    ) 


ii 


=   1  +  0(h) 
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Proof     We  will  use  the  hypotheses  of  the  theorem  to  decompose  S 
into  the  form 

S  =  S  +  h  S  (U5) 

n        n  n 

where  S  is  the  form  of  S  when  constant  steps  are  used  and  S   can  be 

n  n 

bounded.   (This  step  will  be  delayed  until  Lemma  3  below.)   Then  we 
can  examine  I  Is  II  as  follows.   For  fixed  n,  set 

x=Sx  m>n  (U6) 

-m         n  — n  — 

Then 

x   ...    =  S     x 
-m+1         m   -m 

=  Sx+hSx  (1+7) 

-m         m     m   -m 

If  the  constant  step  method  satisfies  the  stability  condition,  then 

Is   "l|  <  krt  .   Lemma  3  will  show  that   Is    <  k, .   Therefore,  Lemma  2 
ii     ii_q  '  '  nM  — •  1 

can  be  applied  to  (^7)  with  k_  =  0 .   Hence, 

k  k  ,T 

M   M    I iom    ll^i     01   II   II 

x    =S   x    <k^e        x 
i  i  — m  i  i    M  n  -n"  -  0        ''-n11 

Hence, 

k  k   T 
l|sm||     :kne01 

I  I    nl  I    _    Q 

so  the  method  satisfies  the  stability  condition. 
Q.E.D. 

Lemma  3   If  a  step  selection  scheme  satisfies  the  hypothesis  of 

Theorem  7,  S   can  be  written  as  S  +  h  S  where  S  is  independent  of 
n  n  n  * 

h  and  n,  and  j  |s    is  uniformly  bounded. 
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Proof 

Part  I :    Interpolation  Technique 
From  (9)  and  (23) 


T  M 
S   =  (I  -  SL      e_)   AC 

n        -n  -2      n 


where  SL      is  independent  of  n.   Now 


Cn  =  diag(l,  hn/hn_15  (hn/hn^)29...) 

=I+diag(0,  \\_±-l,    (h^h^)2-!,...) 

=  I  +  0(h) 
=  I  +  0(h  ) 


since  h   /h         =  1  +  0(h)    and  h  <  h     A-1 
ii     n— jl  —     n 


Therefore , 


Hence, 


S     =    (I    -  I    e^)M  A(I   +   0(h    ))    =  S[I   +   0(h    )] 
n  2  n  n 


S     =  S   0(h    )/h 
n  n       n 


is   uniformly  bounded. 

Part   II :      Variable  Step  Technique 

It   is   shown   in  Tu   [6],   p.    52,   that   the  coefficients  a   and 

have  the  form 

a.        =  a.    +  R      (v    ,.  .  .,  v  ) 

l  ,n  l  a.      n  n-k+2 

l 

3n-      n      =     Sn      +     Rft      (V       '•'•'     Vr,     V  +  9  ^ 

i  ,n  i  p.      n  n-K+£ 

where 


h.   -  h.   , 
v      =     x  J-"1 


1  hi-l 
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where  the  R's  are  rational  functions  that  vanish  at  v   =  v   _  =...=  \)  ,  ,_  -  0 

n    n-1        n-k+2 

(when  the  step  sizes  are  equal),  and  where  the  a.  and  6.  are  the  corresponding 

fixed  step  coefficients.   (This  result  follows  directly  from  the  nonsingular 

system  of  equations  that  determine  a  and  £.)   From  the  hypothesis  of 

Theorem  7,  v.  =  0(h) ,  hence  R  =  0(h)  =  0(h  ).   Therefore,  we  can  express  S  , 

which  involves   the  coefficients  a,    3,   a*  and  3*5as   S  +  h     S  ,  where  S  is  the 

n  n 

matrix  obtained  by  setting  the  a.    equal  to  a. ,  etc.,  and  S   includes  the 

l  ,n  l  n 

terms  R  divided  by  h.   These  can  be  uniformly  bounded. 
Q.E.D. 
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PRACTICAL  IMPLICATIONS 


It  can  be  pointed  out  that  a  common  step  selection  scheme 

satisfies  the  hypotheses  of  Theorem  7 •   If  the  step  h  is  controlled 

so  that  the  local  error  estimate  is  equal  to  e(or  eh  ),  and  if  the  truncation 

n 

r+1  r+2 

error  estimate  is   <J>(t    )   h  +  0(h        )  where   <b(t)   is   the  principle  error 

r     n       n  n  r  ^  * 


function   for  the  method,  we  have 

<|>(t    )    hr+1   +  0(hr+2)    =   e    (or   eh    ) 
n       n  n  n 

Hence , 

h.     =   (-tt|    J      q  +  0(hn)  where  q  =  r+1   or   r, 
n 

provided  that  <j>(t    )    is  not   zero.      If  <\>(t)    is  bounded  away  from  zero  and 

<J>'(t)  is  bounded  above,  then 

n        n+1 

=  1  +  0(h  ) 
n 

and  h  .  /h    >  A  >  0.   Consequently,  we  can  expect  either  method  to  be 
mm  max  — 

stable  if  the  fixed  step  method  is  stable. 

However,  in  practice  we  do  not  send  h  to  the  limit  and  we  often 

do  not  bother  to  adjust  the  step  if  the  change  is  small.   Therefore,  there 

is  reason  to  prefer  the  variable  step  approach.   Brayton,  et.al.  [l]  point 

out  that  the  amount  of  work  in  evaluating  A  y  is  less  in  the  variable  step 

approach  than  in  the  interpolation/Nordsieck  approach,  but  that  the 

coefficients  in  A  and  i      must  be  recomputed  for  each  step.   For  a  large 
n     -n 

system  of  equations,  this  is  not  important  as  A  is  the  same  for  each 
equation  in  the  system.   For  a  few  equations  this  is  a  serious  consideration, 
which  leads  us  to  recommend  variable  step  methods  for  large  systems  of 
equations  and  interpolation  methods  with  the  Nordsieck  vector  for  small  systems 
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